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Abstract 


We  consider  an  alternative  approach  to  the  use  of  nonlinear  stochastic  Markov  processes 
(which  have  a  Fokker-Planck  or  Forward  Kolmogorov  representation  for  density)  in  modeling 
uncertainty  in  populations.  These  alternate  formulations,  which  involve  imposing  probabilis¬ 
tic  structures  on  a  family  of  deterministic  dynamical  systems,  are  shown  to  yield  pointwise 
equivalent  population  densities.  Moreover,  these  alternate  formulations  lead  to  fast  efficient 
calculations  in  inverse  problems  as  well  as  in  forward  simulations.  Here  we  derive  a  class  of 
stochastic  formulations  for  which  such  an  alternate  representation  is  readily  found. 
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1  Introduction  and  Motivation 


In  this  paper  we  consider  classes  of  nonlinear  stochastic  differential  eqnations  (SDE)  with 
corresponding  stochastic  process  density  described  by  Fokker-Planck  (FP)  or  Forward  Kol¬ 
mogorov  (FK)  eqnations  with  nonlinear  drift  (convective  velocity  or  average  transition  rate) 
terms  which  dominate  the  (possibly  nonlinear)  diffusion  (variability  in  transition  rate)  terms. 
We  develop  a  general  algorithmic  approach  for  converting  these  computationally  difficult 
nonlinear  SDF  to  an  equivalent  (in  a  sense  to  be  made  precise  below)  probabilistic  formu¬ 
lation  which  is  much  more  amenable  to  fast  (and  parallel)  computations.  A  major  feature 
of  our  approach  here  is  the  bidirectional  nature  of  our  derivations.  That  is,  we  show  how 
to  transform  from  a  given  SDF  to  the  corresponding  equivalent  probabilistic  formulation, 
and  from  a  given  probabilistic  formulation  to  the  corresponding  SDE.  Our  results  are  pre¬ 
sented  in  two  distinct  steps.  First,  we  show  an  equivalence  for  several  classes  of  nonau- 
tonomous  affine  differential  equations  (in  both  scalar  and  system  forms).  Then  we  extend 
theses  equivalences  to  rather  general  classes  of  nonlinear  differential  equations  using  invert¬ 
ible  transformation  techniques  between  nonlinear  differential  equations  and  the  class  of  affine 
differential  equations  introduced  in  the  first  step.  We  illustrate  the  ideas  with  several  dif¬ 
ferent  nonlinear  examples  including  two  important  examples  with  growth  or  transition  rates 
encountered  frequently  in  modeling  populations,  label  decay,  tumor  growth,  etc.  These  are 
the  popular  Verhulst-Pearl  logistic  growth  rate  g{x)  =  rx  {l - ^  and  the  general  tran¬ 

sition  rates  g{x,t)  =  (ao(t)  —  ai{t)\\ix)x  of  which  the  standard  Gompertz  growth  rates 
g{x)  =  rln  ^  x  =  r(lnfi;  —  lna;)a;  are  special  cases.  Motivating  ideas  and  a  brief  recap  of 

results  from  previous  efforts  on  population  models  are  first  summarized  here  and  in  the  next 
section. 

A  fundamental  modeling  construction  in  many  areas  of  science  is  the  nonlinear  Markov 
process  as  characterized  by  discrete  or  continuous  time  with  discrete  (e.g.,  Poisson)  or  con¬ 
tinuous  (e.g.,  diffusion)  state  processes  [2,  22,  24,  27,  31,  33].  Here  we  focus  on  modeling  with 
a  general  nonlinear  Markov  diffusion  process  with  finite  mean  and  variance.  Mathematically 
this  leads  to  a  stochastic  differential  equation  (SDE)  of  the  form 

dX{t)  =  g{X{t),t)dt  +  a{X{t),t)dW{t), 

which  in  turn  can  be  investigated  with  Fokker-Planck  (FP)  or  Forward  Kolmogorov  equations 

f  f)  1 

^(a^G)  +  ■^{9{x,t)u{x,t))  =  -  —  {a^{x,t)u{x,t)),  (1.1) 

for  the  corresponding  probability  density  u.  These  FP  models  are  ubiquitous  in  mathematics 
and  physics  (e.g.,  particle  transport,  hltering),  biology  (population  models),  hnance  (e.g., 
Black-Scholes  equations)  among  other  areas  of  applications  [2,  22,  27,  29].  In  many  of  these 
applications  one  has  what  is  commonly  referred  to  as  convection  dominated  diffusion  which 
occurs  when  we  have  g  »  a.  In  this  case  the  Fokker-Planck  equations  are  notoriously 
difficult  to  solve  especially  when  g  depends  on  time.  This  is  a  serious  drawback  in  forward 
simulations  but  can  be  untenable  in  inverse  problem  calculations  [20]. 
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In  these  applications  one  begins  with  a  stochastic  Markov  diffusion  process  X{t)  with  real¬ 
izations  x(t)  representing  a  structure  variable  (level)  such  as  size  (length,  weight,  volume, 
etc.),  label  intensity,  chronological  or  physiological  age,  spatial  location,  etc.,  that  changes 
according  to  a  mean  rate  law  g{x,t)  with  mean  variance  a‘^{x,t).  Such  models  also  arise 
in  population  biology  with  class  structure  [4,  5,  6,  7,  9,  10,  12,  14],  complex  nodal  network 
models  (in  network  security,  social/insurgency  networks,  logistic  and  production  networks) 
[3,  16],  fluorescence  intensity  of  labeled  proliferating  cell  populations  [17,  28]  and  general 
hyperbolic  transport  systems  in  random  or  uncertain  environments. 

To  incorporate  uncertainty  or  variability  into  structured  deterministic  dynamical  models, 
several  approaches  have  been  considered  in  the  literature.  Of  those  of  interest  to  us  here,  one 
involves  a  stationary  probabilistic  structure  on  a  family  of  structured  deterministic  dynamical 
systems,  while  the  other  is  constructed  based  on  the  assumption  that  movement  from  one 
structure  level  to  another  can  be  described  by  a  stochastic  diffusion  process.  As  noted  above, 
for  computational  purposes  the  latter  can  be  represented  mathematically  by  a  Fokker-Planck 
equation.  Even  though  these  two  formulations  are  conceptually  quite  different,  in  this  paper 
we  show  for  a  class  of  examples  (both  linear  and  nonlinear  differential  equations)  that  they 
are  equivalent  in  terms  of  corresponding  probability  density  functions.  Numerical  methods 
for  the  probabilistic  formulation  are  quite  fast  and  lead  to  alternative  methods  for  solution  of 
the  Fokker-Planck  or  Forward  Kolmogorov  equation  associated  with  the  stochastic  diffusion 
process.  Thus  we  establish  that  there  are  several  classes  of  Fokker-Planck  inverse  problems 
(which  are  computationally  intensive)  that  can  be  readily  converted  to  inverse  problems  for 
probabilistic  structures  on  deterministic  systems  which  can  be  solved  efficiently. 

Because  we  are  primarily  interested  in  modeling  transition  uncertainty  in  this  paper,  for 
simplicity,  we  will  not  consider  either  sink  or  source  terms  in  our  formulations. 

Unless  otherwise  indicated,  a  capital  letter  is  used  to  denote  a  random  variable  throughout 
the  presentation,  and  a  corresponding  small  letter  to  denote  its  realization.  We  use  A/’(/i,  cr^) 
to  denote  a  normal  distribution  with  mean  p  and  variance  A/'(/x,  S)  for  a  multivariate 
normal  distribution  with  mean  vector  /x  and  covariance  matrix  S,  and  E(Z)  for  the  expecta¬ 
tion  of  random  variable  Z.  We  next  give  careful  detailed  formulations  of  the  two  approaches 
we  consider. 


1.1  Stochastic  Formulation 

A  stochastic  formulation  can  be  motivated  by  recognizing  that  local  factors  (such  as  environ¬ 
mental  or  emotional  fluctuations)  can  have  a  signihcant  influence  on  the  individual  transition 
or  transfer  rates  from  one  structure  level  to  another.  For  example,  in  [5,  7,  9,  10,  20]  the 
growth  rate  of  two  different  marine  species  (mosquitofish  and  shrimp)  are  affected  by  sev¬ 
eral  environmental  factors  such  as  temperature,  dissolved  oxygen  level  and  salinity.  In  such 
examples,  the  stochastic  formulation  is  constructed  under  the  assumption  that  movement 
from  one  structure  level  (size  in  these  examples)  to  another  can  be  described  by  a  stochastic 
diffusion  process  [2,  20,  24,  31].  Let  {X(t)  :  t  >  0}  be  a  Markov  diffusion  process  with  X(t) 
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representing  structure  level  at  time  t  (i.e.,  each  process  realization  corresponds  to  the  struc¬ 
ture  trajectory  of  an  individual).  Then  X{t)  is  described  by  the  Ito  stochastic  differential 
equation  (we  refer  to  this  equation  as  the  stochastic  rate  model  (SRM)) 

dX{t)  =  g{X{t),t)dt  +  a{X{t),t)dW{t),  (1.2) 

where  W (t)  is  the  standard  Wiener  process  [2,  24].  Here  g{x,  t)  denotes  the  average  transition 
rate  (the  hrst  moment  of  rate  of  change  in  structure  level)  of  individuals  with  structure  level 
X  at  time  t,  and  is  given  by 

g{x,t)  =  ^lim^^E{AX(t)|X(f)  =  x}  .  (1.3) 

The  function  a{x,  t)  represents  the  variability  in  the  transition  rate  of  individuals  (the  second 
moment  of  rate  of  change  in  structure  level)  and  is  given  by 

(r^{x,t)  =  Jim^^E{[AX(t)]2|X(t)  =  x}  .  (1.4) 

Hence,  the  transition  process  (growth  in  size  in  the  marine  examples  cited  above)  of  each 
individual  is  stochastic,  and  each  individual  changes  structure  level  according  to  the  stochas¬ 
tic  rate  model  (1.2).  Thus,  for  this  formulation  the  transition  uncertainty  is  introduced  into 
the  entire  population  by  the  stochastic  transition  of  each  individual.  In  addition,  individuals 
with  the  same  structure  level  at  the  same  time  have  the  same  uncertainty  in  transition, 
and  individuals  also  have  the  possibility  of  reducing  their  structure  level  during  a  transition 
period. 

With  this  assumption  on  the  transition  process,  we  obtain  the  Fokker-Planck  (FP)  or  forward 
Kolmogorov  model  for  the  population  density  u,  which  was  carefully  derived  in  [31]  among 
numerous  other  places  and  subsequently  studied  in  many  references  (e.g.,  [2,  20,  24]).  The 
equation  and  appropriate  boundary  conditions  are  given  by 

du  d  1  (9^ 

-^{x,t)  + -^{g{x,t)u{x,t))  = -  —  {a^{x,t)u{x,t)),  xE{0,L),  t>0, 

g{0,t)u{0,t)  -  ^-^{a^{x,t)u{x,t))\^=o  =  0, 

1  d 

g{L,t)u{L,t)  -  -  —  {a^{x,t)u{x,t))\^=L  =  0, 
u{x,  0)  =  Uo{x). 

Here  L  is  the  maximum  structure  level  that  individuals  may  attain  in  any  given  time  period. 
Observe  that  the  boundary  conditions  in  (1.5)  provide  a  conservation  law  for  the  FP  model. 
Because  both  sink  (death)  and  reproduction  (birth)  rates  are  assumed  zero,  the  total  number 

.L 

in  the  population  is  a  constant  given  by  /  uo{x)dx.  In  addition,  we  observe  that  with  the 

Jo 

zero-flux  boundary  condition  at  zero  (minimum  structure  level)  one  can  equivalently  set 
X{t)  =0  if  X{t)  <0  for  the  stochastic  transition  model  (1.2)  in  the  sense  that  both  are  used 
to  keep  individuals  in  the  system  at  the  boundary  levels.  This  means  that  if  the  structure 
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level  of  an  individual  is  decreased  to  the  minimum  value,  it  remains  in  the  system  at  that 
level  with  the  possibility  to  once  again  increase  its  level. 

As  we  have  already  mentioned,  the  above  equations  present  formidable  computational  chal¬ 
lenges  when  drift  (or  convection)  g  dominates.  While  this  is  certainly  true  for  forward 
simulations,  it  is  even  more  of  a  challenge  in  the  inverse  problems  considered  in  [20].  Thus 
it  is  of  great  interest  to  explore  alternative  formulations  that  might  lead  to  modeling  repre¬ 
sentations  which  represent  essentially  the  same  transport  process  at  the  density/distribution 
level  and  that  involve  a  much  more  efficient  computational  implementation.  To  this  end  in 
the  next  section  we  describe  an  alternative,  but  in  a  certain  sense  (to  be  made  more  precise 
below)  equivalent,  formulation  for  such  dynamic  processes  with  uncertainty  in  transition. 


1.2  An  Alternative  Probabilistic  Formulation 

To  offer  an  alternative  formulation  to  stochastic  Markov  processes  for  populations  with 
uncertainty  in  transition  rates,  we  turn  to  population  models  (labeled  proliferating  cells  [17], 
mosquitohsh  [7,  14],  shrimp  [5,  9,  10])  where  we  have  used  structured  dynamical  models  with 
uncertainty.  The  probabilistic  formulation  we  present  is  motivated  by  the  observation  that 
in  populations,  intrinsic  variables  (genetic  differences  or  non-lethal  infections  of  some  chronic 
disease)  can  have  an  effect  on  individual  class  transition  rates.  For  example,  in  many  marine 
species  such  as  mosquitohsh,  females  grow  faster  than  males,  which  means  that  individuals 
with  the  same  size  may  have  different  growth  rates.  In  labeled  cell  populations  where  label 
intensity  is  the  structure  variable,  label  decay  rates  may  vary  across  individuals  during 
cell  proliferation  due  to  variable  decay  in  cell  cycle  stages.  The  probabilistic  formulation 

is  constructed  based  on  the  assumption  that  each  individual  does  change  according  to  a 

dx 

deterministic  rate  model  as  posited  in  the  Sinko-Streifer  [32]  formulation, 

but  that  different  individuals  may  have  different  structure-dependent  rates.  Based  on  this 
underlying  assumption,  one  partitions  the  entire  population  into  (possibly  a  continuum  of) 
subpopulations  where  individuals  in  each  subpopulation  have  the  same  structure-dependent 
transition  rates,  and  then  assigns  a  probability  distribution  to  this  partition  of  possible 
transition  rates  in  the  population. 

To  be  more  precise  here  we  can  describe  this  construction  in  terms  of  population  growth 
models  where  size  (length  in  mosquitohsh,  weight  in  shrimp)  is  the  structure  variable.  Then 
the  growth  process  for  individuals  in  a  subpopulation  with  rate  g  is  assumed  to  be  described 
by  the  dynamics 

=  9{x{t;g),t),  geG,  (1.6) 

where  ^  is  a  collection  of  admissible  rates.  Model  (1.6)  combined  with  the  probability 
distribution  imposed  on  G  will  be  called  the  prohabilistie  rate  distribution  (PRD)  model  in 
this  paper.  Hence,  we  can  see  that  for  the  probabilistic  formulation,  the  rate  uncertainty  is 
introduced  into  the  entire  population  by  the  variability  of  transition  (growth)  rates  among 
subpopulations.  In  the  literature  for  size-structured  population  models,  it  is  common  to 
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assume  that  growth  rate  is  a  nonnegative  function,  that  is,  no  loss  in  size  occurs.  However, 
individuals  may  experience  loss  in  size  due  to  disease  or  some  other  involuntary  factors. 
Hence,  we  will  permit  these  situations  in  this  formulation,  but  for  simplicity  we  assume  that 
growth  rate  in  each  subpopulation  is  either  a  nonnegative  function  or  a  negative  function, 
that  is,  the  size  of  each  individual  may  be  either  nondecreasing  or  decreasing  continuously 
in  its  growth  period. 

With  this  assumption  of  a  family  of  admissible  transition  rates  and  an  associated  probability 
distribution  (i.e.,  the  PRD  model),  one  can  obtain  a  generalization  of  the  Sinko-Streifer  model 
for  densities  as  a  function  of  time  and  structure  variables.  (This  was  called  the  growth  rate 
distribution  or  (GRD)  model  as  formulated  and  studied  in  [4,  6,  7,  12,  14].)  The  model, 
which  here  will  be  called  the  class  rate  distribution- Sinko-Streifer  (CRDSS)  model,  consists 
of  solving 

vt{x,t;g)  +  {g{x,t)v{x,t;g))^  =  0,  a;e(0,  L),  t  >  0, 

g{0,  t)v{0,  t;g)  =  0iig>0  or  g{L,  t)v{L,  t;g)  =  0ii  g  <0,  (1.7) 

v{x,0-,g)  =  vo{x-,g), 

for  a  given  g  E  Q  and  then  “summing”  (with  respect  to  the  probability)  the  corresponding 
solutions  over  all  g  E  Q.  If  v{x,t-,g)  is  the  population  density  of  individuals  with  class 
structure  value  x  at  time  t  having  transition  rate  g,  the  expectation  of  the  total  population 
density  for  class  x  at  time  t  is  given  by 

u{x,t)=[  v{x,t;g)dV{g),  (1.8) 

Jgeg 

where  P  is  a  probability  measure  on  Q.  Hence  the  CRDSS  model  consists  of  the  Sinko- 
Streifer  equation  (1.7)  with  (1.8)  to  compute  the  population  density  u{x,t).  Thus,  this 
probabilistic  formulation  involves  a  stationary  probabilistic  structure  on  a  family  of  deter¬ 
ministic  dynamical  systems,  and  V  is  the  fundamental  “parameter”  that  is  to  be  estimated 
by  either  parametric  or  nonparametric  methods  (which  depends  on  the  prior  information 
known  about  the  form  for  V).  As  detailed  in  [7,  12],  this  class  rate  distribution  model  is 
sufficiently  rich  to  exhibit  a  number  of  phenomena  of  interest,  for  example,  dispersion  and 
development  of  two  modes  from  one.  Moreover,  of  paramount  importance  to  us  here,  this 
formulation  offers  tremendous  computational  advantages  in  that  it  is  what  may  be  accurately 
termed  embarrassingly  parallel  [14]. 

Observe  that  if  all  the  subpopulations  have  nonnegative  rate  functions,  then  we  need  to  set 
g{L,t)v{L,t]  g)  =  0  for  each  g  E  Q  m.  order  to  provide  a  conservation  law  for  the  CRDSS 
model.  Specihcally  if  L  denotes  the  maximnm  attainable  class  valne  of  individnals  in  a  life 
time,  then  it  is  reasonable  to  set  g{L,  f)  =  0  (as  commonly  done  in  the  literature).  However,  if 
we  just  consider  the  model  in  a  short  time  period,  then  we  may  choose  L  snfficiently  large  so 
that  v{L,  t)  is  negligible  or  zero  if  possible.  We  observe  that  if  there  exist  some  snbpopnlations 
whose  rates  are  negative,  then  we  can  not  provide  a  conservation  law  for  these  snbpopnlations 
as  5f(0,  t)  <  0.  Hence,  in  this  case,  once  the  class  valne  of  an  individnal  is  decreased  to  below 
the  minimnm  valne,  then  that  individnal  will  be  removed  from  the  system.  In  other  words. 
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we  exclude  those  individuals  whose  class  value  go  below  the  minimum  size.  This  effectively 
serves  as  a  sink  for  these  subpopulations. 


2  Summary  of  Previous  Findings 


We  compare  here  the  probabilistic  formulation  approach  to  incorporating  the  class  rate 
uncertainty  into  a  structured  population  model  with  the  stochastic  formulation  of  Section 


1.1. 


The  discussions  in  Sections  1.1  and  1.2  indicate  that  these  two  formulations  are  conceptually 
quite  different.  One  entails  imposing  a  probabilistic  structure  on  the  set  of  possible  transition 
rates  permissible  in  the  entire  population  while  the  other  involves  formulating  transition  as 
a  stochastic  diffusion  process.  However,  the  analysis  in  [9]  reveals  that  in  some  cases  the 
structure  distribution  (the  probability  density  function  of  X{t))  obtained  from  the  stochastic 
rate  model  is  exactly  the  same  as  that  obtained  from  the  PRD  model.  For  example,  if  we 
consider  the  two  models 


stochastic  formulation:  dX(t)  =  bo{X(t)  +  co)dt  +  V^ao{X(t)  +  co)dW(t) 

dx(t;  h) 


probabilistic  formulation: 


=  (6  -  aot)(a;(t;6)  +  Co), 


dt 

B~U{ba,al), 


(2.1) 


and  assume  their  initial  structure  distributions  are  the  same,  then  we  obtain  at  each  time 
t  the  same  structure  distribution  from  these  two  distinct  formulations.  Here  <^o  and  Cq 
are  positive  constants  (for  application  purposes),  and  R  is  a  normal  random  variable  with 
b  a  realization  of  B.  Moreover,  by  using  the  same  analysis  as  in  [9]  we  can  show  that  if  we 
compare 


stochastic  formulation:  dX{t)  =  {bo  +  (Tgt)(X(t)  +  co)dt  +  \/^ao{X{t)  +  Co)dW{t) 


probabilistic  formulation: 


dx{t]  b) 
dt 


b{x{t;  b)  +  Co),  6  G  M  with  B  ~  M{bo,  Ug) 


(2.2) 


with  the  same  initial  structure  distributions,  then  we  can  also  obtain  at  each  time  t  the 
same  structure  distribution  for  these  two  formulations.  In  addition,  we  see  that  both  the 
stochastic  rate  models  and  the  probabilistic  rate  models  in  (2.1)  and  (2.2)  reduce  to  the  same 
deterministic  growth  model  x  =  bo{x  +  co)  when  there  is  no  uncertainty  or  variability  in  rate 
(i.e.,  (To  =  0)  even  though  both  models  in  (2.2)  do  not  satisfy  the  mean  rate  dynamics 


=  (,„(E(2f  (())  +  c„)  (2.3) 

while  both  models  in  (2.1)  do.  This  last  observation  was  critical  in  the  early  efforts  of  [9,  10] 
which  were  derived  under  the  additional  constraint  that  (2.3)  must  hold.  This  was  motivated 
by  available  shrimp  data  of  longitudinal  measurements  of  average  shrimp  weight  (in  gms). 
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i.e.,  an  observation  of  x{t)  =  E(X(t)).  In  this  earlier  work  it  was  found  that  an  affine  growth 

law  =  g{x{t))  =  bo{x{t)  +  Cq)  yields  a  good  £t  to  this  data  for  early  shrimp  growth. 

This  led  to  a  search  for  equivalent  mathematical  representations  which  also  satisfied  this 
extra  condition. 

More  specifically,  one  can  prove  that  the  formulations  in  (2.1)  generate  stochastic  processes 
X{t)  which  both  satisfy  the  mean  rate  dynamics  (2.3)  and  yield  processes 

X{t)  =  -co  +  {Xo  +  co)Y{t), 

where 

YpRoit)  =  exp(5t  -  where  B  ~  A/'(6o,  (^l),  (2.4) 

XsRM  it)  =  exp  ^{bot  -  +  *^0  y  y/^dW{T)^  .  (2.5) 

Moreover  it  was  shown  that  for  each  time  t,  both  YpR£,{t)  and  Yspuit)  are  normally  dis¬ 
tributed  with  identical  means  and  variances.  Thus  under  the  additional  reasonable  assump¬ 
tion  (trivially  true  for  non-random  initial  condition)  that  the  random  variables  Xq  and  each 
of  YpRpit)  and  YsRMit)  are  independent  we  find  that  each  of  the  stochastic  processes  de¬ 
rived  from  (2.1)  possess  at  each  time  t  the  same  distribution.  That  is,  at  each  time  t  each  of 
the  processes  X(t)  have  the  same  probability  density.  Finally,  the  two  stochastic  processes 
are  NOT  the  same.  This  can  be  seen  immediately  from  (2.4)  and  (2.5),  but  also  from  a 
direct  calculation  of  the  covariances  for  Yppp  and  Ysrm  which  we  shall  carry  out  below.  In 
summary,  while  the  two  formulations  of  (2.1)  generally  lead  to  different  processes,  one  can 
argue  that  they  are  equivalent  in  the  sense  that  they  possess  the  same  probability  density 
at  any  time  t.  For  the  subsequent  discussions  in  this  presentation,  we  shall  refer  to  this  as 
pointwise  equivalence  in  density.  This  density  must  satisfy  the  corresponding  Fokker-Planck 
or  Forward  Kolmogorov  equation  for  the  stochastic  formulation  in  (2.1).  Thus  if  one  wishes 
to  obtain  a  numerical  solution  of  such  a  Fokker-Planck  equation,  one  possibility  is  to  consider 
the  numerical  solution  of  the  equivalent  but  more  readily  solved  CRDSS  formulation  of  (2.1). 
For  the  particular  forms  of  (2.1)  and  (2.2),  this  approach  was  demonstrated  to  be  a  compu¬ 
tationally  advantageous  strategy  in  [11].  These  findings  lead  to  a  natural  research  question: 
Are  there  general  classes  of  Fokker-Planck  equations  that  can  be  converted  to  an  equivalent 
(in  the  distributional  sense  described  above)  CRDSS  which  can  be  efficiently  solved  numer¬ 
ically  for  the  desired  probability  density  function?  A  positive  answer  to  this  question  is  the 
primary  focus  of  this  paper  and  results  are  given  in  the  next  sections.  In  particular  in  the 
next  section  we  develop  general  techniques  to  show  equivalence  for  large  classes  of  affine 
differential  equations. 


3  Equivalence  between  Probabilistic  and  Stochastic  For¬ 
mulations  with  Affine  Dynamics 


In  this  section,  we  turn  to  several  class  of  examples  with  affine  dynamics  for  which  one 
can  establish  the  desired  equivalence  between  the  probabilistic  and  stochastic  formulations 
given  above.  The  probabilistic  formulations  we  consider  here  involve  a  hnite-dimensional 
parameter  family  of  structure  rates  of  change.  For  example,  if  the  probabilistic  formulation 
is  governed  by  scalar  differential  equations,  then  we  assume  that  all  the  subsystems  have  the 
same  functional  form  g{x,  t]  b)  for  the  structure  rates  of  change  but  the  values  of  parameters 
b  =  {bo,  bi, . . . ,  6n-i)^  vary  across  the  system. 

In  establishing  our  results  (and  to  discuss  corresponding  covariances  below),  the  following 
relationship  between  normal  distribution  and  log-normal  distribution  [21,  page  109]  is  heavily 
used. 


Lemma  3.1  IflnZ  ~  J\f  {jj,,  a'^),  then  Z  is  log-normally  distributed,  where  its  probability 
density  function  fz{z)  is  defined  by 


fz{z) 


zy/^a 


exp 


V  )' 


and  its  mean  and  variance  are  given  as  follows 


E(Z)  =  exp(/i  -|-  Var{Z)  =  [exp((T^)  —  1]  exp(2/i  -|-  a^). 


In  our  subsequent  arguments  we  shall  also  need  the  following  basic  result  on  the  process 
generated  by  Ito  integrals  of  Wiener  processes  that  can  be  found  in  [27,  Sec  4.3,  Thm  4.11]. 


Lemma  3.2  Let  T  be  any^  positive  constant.  Then  for  a  non-random  function  f  G  L^{0,  T), 
the  Ito  integrals  Q{f)  =  /  f{s)dW (s)  for  0  <  t  <T  yield  a  Gaussian  stochastic  process  with 

Jo  ^  ^ 

pointwise  distributions  Af  J  f^{s)ds^.  Moreover,  Cov{Q{t),Q{t-\-T))  =  j  p{s)ds  for 
all  r  >  0. 


3.1  Scalar  Differential  Equations  -  Case  I 

In  the  hrst  case  we  derive  conditions  under  which  the  probabilistic  and  stochastic  formula¬ 
tions  generate  stochastic  processes  with  the  same  distributions  (normal  in  the  case  the  initial 
condition  is  a  fixed  constant)  at  each  time  t. 
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Probabilistic  formulation  The  probabilistic  formulation  considered  has  the  following 
form 


dx(t;  b) 
dt 


=  a{t)x{t]  b)  +  7(t)  +  b  ■  Q{t),  b  =  (6o,  &i,  •  •  • ,  hn-iY  G 


Bj  ~  CTj),  j  =  0, 1,  2, . . . ,  n  —  1,  which  are  mutually  independent, 


(3.1) 


where  a,  7  and  q  =  (po,  Pi,  ■  ■  ■ ,  Pn-i)^  are  non-random  functions  of  t,  and  b  is  chosen  as 
realizations  of  B  =  (Bq,  i?i, . . . ,  Bn-i)'^ ■  Hence,  the  dynamics  of  an  individual  with  initial 
condition  xq  in  a  subsystem  with  its  rates  of  change  having  parameter  values  b  is  described 
by  the  deterministic  model 


dx(t]  b) 
dt 


=  a{t)x{t;  b)  -|-  7(t)  -|-  b  ■  Q{t),  a;(0)  =  Xq. 


Multiplying  both  sides  of  the  above  equation  by  the  integrating  factor  exp 
we  hnd  that 


a{s)ds 


a; exp  (  —  /  a{s)ds 


=  7(t)  exp  (  —  /  a{s)ds  )  -|-  b  ■  git)  exp  (  —  /  a{s)ds 


d 
dt 

We  then  integrate  both  sides  of  the  above  equation  to  obtain 


x{t]XQ,h)  =  xoexpj  /  Q;(s)(is  )  -|-  /  7(5)  exp 

Jo 

r't 


a{T)dT  ]  ds 


+b 


Q[s) exp 


a{T)dT  )  ds. 


(3.2) 


We  assume  that  all  the  subsystems  have  the  same  probability  density  function  for  initial 
condition  Xq,  which  is  independent  of  Bj,  j  =  0, 1,  2, . . . ,  n  —  1.  Let  X(t)  =  xit]  Xq,  B)  and 


Y{t)=  /  7(s)exp 
Jo 

Then  we  have  that 


a{T)dT  )  ds  -|-  B  •  /  g{s)  exp  (  /  a{T)dT  ]  ds. 


X(t)  =  Xoexp  (  /  Q;(s)ds  )  -|- F(t) 


(3.3) 

Note  that  Bj  ~  A/’(/ij,  (j|),  and  Bj,  j  =  0, 1,  2, . . . ,  n  —  1,  are  mutually  independent.  Hence, 
we  hnd  that  for  any  hxed  t,  Y (t)  is  normally  distributed  with  mean  dehned  by 


(7(5)  -|-  n  ■  q{s))  exp  (  /  a{T)dT  ]  ds. 


(3.4) 


where  fi  =  {/xq,  fxi, . . . ,  jXn-iY' ,  and  variance  dehned  by 


/  f?j(s)exp(  /  a{T)dT\  ds 


(3.5) 
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Stochastic  formulation  Based  on  the  above  discussions  of  the  probabilistic  formulation 
(3.1)  we  see  that  if  all  the  individuals  in  the  entire  system  have  the  same  fixed  initial  condition 
Xq,  then  X{t)  is  also  normally  distributed  for  any  fixed  time  t.  Based  on  this  piece  of 
information,  the  stochastic  model  is  chosen  to  have  the  form 


dX{t)  =  [a{t)X{t)  +  ^{t)]dt  +  7]{t)dW{t),  X(0)  =  Xq,  (3.6) 

where  a,  ^  and  77  are  non-random  functions  of  t,  and  Xq  is  independent  of  W {t).  Multiplying 


both  sides  of  the  above  equation  by  exp 


a{s)ds  ]  we  hnd  that 


X  (t)  exp 


a{s)ds 


=  ^(^)exp 


a{s)ds  )  dt  +  rj{t)  exp 


Integrating  both  sides  of  the  above  equation  we  obtain  that 


X(t)=Xoexp(  /  Q;(s)(is  ]  -|- F(t), 


a{s)ds^  dW(t). 


(3.7) 


where  Y  (t)  is  dehned  by 


a 


Y{t)=  /  e(5)exp 
Jo 

By  Lemma  3.2,  for  any  hxed  t  we  have  that 

J  ri{s)exp^—J  Q;(r)(ir^  (ihL(s)  ~  A/”  (  0, 

Thus,  we  hnd  that 

Y{t)  ~  Af  (f  ^(s)exp  (f  a{T)dT  ]  ds. 


(r)dr^  ds-fexp  a{T)dT^  J  r]{s)exp^—J  Q;(r)dr^  dhL(s) 


rt  r 


77(5)  exp  (  —  /  a{T)dT 


ds 


rt  V 


r]{s)  exp  (  /  a{T)dT 


n  2 


ds 


(3.8) 


Equivalent  conditions  By  (3.3),  (3.4),  (3.5),  (3.7)  and  (3.8),  we  know  that  if  functions 
77  and  Qj,  and  constants  /ij,  aj  and  n  satisfy  the  following  two  equalities 

J  ^{s)  exp  a{T)dT^  ds  =  j  [^{s)  + /j,  ■  Q{s)]exp  a{T)dT^  ds  (3.9) 


and 


/t  ~  /  \'^  Til  ~ 

77(5)  exp  (  /  a{T)dT\  ds  =  '^^a‘^  /  f7j(s)exp(  /  a{T)d 


T  1  ds 


(3.10) 


then  the  probabilistic  formulation  (3.1)  and  the  stochastic  formulation  (3.6)  yield  stochastic 
processes  that  are  pointwise  equivalent  in  density. 


Based  on  the  equivalent  conditions  (3.9)  and  (3.10),  we  will  derive  the  specihc  forms  of  the 
functional  parameters  and/or  values  of  the  parameters  of  the  corresponding  pointwise  equiv¬ 
alent  (stochastic/probabilistic)  formulation  in  terms  of  the  functional  parameters  and/or 
parameters  of  the  known  (probabilistic/stochastic)  formulation. 
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3.1.1  Probabilistic  Formulation  to  Stochastic  Formulation 


Here  we  assume  that  probabilistic  formulation  (3.1)  is  known,  and  we  want  to  determine  its 
corresponding  stochastic  formulation.  In  other  words,  we  need  to  determine  functions  ^  and 
f]  in  terms  of  functions  7  and  gj,  and  constants  fij,  Cj  and  n.  By  (3.9),  it  is  obvious  that  if 
function  ^  is  chosen  to  be 

^(t)  =^(t)  +  /X-  Q{t), 

then  (3.9)  holds.  Differentiating  both  sides  of  (3.10)  with  respect  to  t  yields  that 


n—1 


v\t)  +  2a{t)J2 


j=0 


gj{s)exp  (  /  a{T)dT  )  ds 


n—1 


2  Ed 

i=o 


gj{t)  +  a{t)  /  f?j(s)exp(  /  Q;(r)(ir  )  ds 


'0 


f7j(s)exp  (  /  a{T)dT  )  ds. 


which  can  be  simplihed  as  follows 

n—1 


,m=2j2a]  gj(t)  j  gj{s)exp^J  Q;(r)dr^  ds. 

Hence,  if  q  is  assumed  to  have  the  property  that  the  right-hand-side  of  the  above  equation 
is  nonnegative  for  any  t  >  0,  then  we  can  choose  7  to  be 


r]{t)  = 


pt  /  pt 

2  y  J 


a{T)dT  )  ds 


so  that  (3.10)  holds. 

3.1.2  Stochastic  Formulation  to  Probabilistic  Formulation 

Next  we  assume  that  stochastic  formulation  (3.6)  is  known,  and  we  wish  to  determine  its 
corresponding  probabilistic  formulation.  In  other  words,  we  need  to  determine  functions  7 
and  pj,  and  constants  pj,  aj  and  n  in  terms  of  fnnctions  ^  and  p.  It  is  obvious  that  if  we  set 

7(f)  +  p  ■  e(t)  =  ^(t),  (3.11) 

then  (3.9)  holds.  By  (3.11)  and  (3.10)  we  know  that  we  have  numerous  different  choices  for 
the  probabilistic  formnlation.  Here  we  illustrate  two  of  them. 

Choice  1  One  of  the  simple  choices  is  to  choose 

n  =  2,  7  =  0,  po  =  1,  Pi  =  0,  go{t)=^{t).  (3.12) 
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Then  by  the  above  equalities  and  (3.10)  we  have 


/  ^?l(’S)exp  (  /  a{T)dT  )  ds 


=  /  ri^{s)exp  I  2  /  a{T)dT  \  ds  —  aQ  /  .^(s)exp  I  /  a{T)dT  |  ds  , 
Jo  \  Js  J  iJo  \Js  J  \ 


(3.13) 


which  implies  that  we  need  to  choose  ao  sufficiently  small  such  that  its  right-hand  side  is 
greater  than  0.  Now  by  (3.13)  we  have 

<^i  [  ^?l('S)exp^  [  Q;(r)(ir^  ds 


rf‘{s)  exp  (  2  /  a{T)dT  \  ds  —  i  ^{s)  exp  (  /  a{T)dT  )  ds 


Differentiating  both  sides  of  the  above  equation  with  respect  to  t  we  obtain  that 


<^l^?l(^)  ~  IT  /  I  2  /  a{T)dT  I  ds  —  (Tq  I  /  .^(s)exp  I  /  a(T)dT  |  ds 

Jo  \  Js  J  \Jo  \Js  J 


-a{t)  /  77^(s)exp  I  2  /  a{T)dT  j  ds  —  cig  |  /  .^(s)exp  |  /  a{T)dT  j  ds 

Jo  \  Js  J  \Jo  \Js  J 


Hence,  we  can  just  assign  any  positive  value  for  ai,  and  then  use  the  above  equality  to 
determine  function  qi.  Hence,  one  simple  choice  of  the  corresponding  pointwise  equivalent 
probabilistic  formulation  for  (3.6)  is  as  follows 

^  a(t)a;(t;  6o,  &i)  +  &o^(t)  +  feiPi(t),  6o,  e  M, 

Ho  ~  ■A/”(l,  (Jq),  Hi  ~  A/”(0,  (T^), 

where  Uo,  Ui  and  pi  are  chosen  according  to  the  above  discussions. 


Choice  2  Another  simple  choice  is  to  choose 

n  =  l,  7(t)=^(t),  /io  =  0,  0-0  =  1. 

Thus,  by  the  above  equalities  and  (3.10)  we  hnd 


I  Po('S)exp  Q;(r)drJ  ds  =  ^y  7^(s)exp  ^2  J  a{T)dT  j  ds  j  . 


(3.14) 
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Differentiating  both  sides  of  the  above  equation  with  respect  to  t  we  have  that 

d 


Qo{t)  = 


dt 


ri‘^{s)exp  (  2  /  a{T)dT  ds 


-a{t) 


Uo 


if  {s)  exp  (  2  /  a{T)dT 


(3.15) 


Hence,  another  simple  choice  of  the  corresponding  pointwise  equivalent  probabilistic  formu¬ 
lation  for  (3.6)  is  as  follows 


dx(t;  h) 


dt 

B  ~  A/'(0, 1), 


=  a{t)x{t]h)  +  ^{t) +  hQQ{t),  beR, 


where  Qq  is  dehned  by  (3.15). 


3.2  Scalar  Differential  Equations  -  Case  II 

In  this  section,  we  consider  examples  where  the  probabilistic  and  stochastic  formulations  lead 
to  processes  that  have  pointwise  equivalent  densities  (either  log-normal  or  shifted  log-normal, 
depending  on  the  value  of  the  parameter  (see  details  below),  when  the  initial  condition  is  a 
hxed  constant). 


Probabilistic  formulation  The  probabilistic  formulation  considered  has  the  form 
dxit]  b) 


dt 


=  (b  ■  p{t)  +  7(t))(a;(t;  b)  +  c),  b  =  {bo,  bi,...,  6„_i)  G 


(3.16) 


B,  M{pj,  =  0,l,2,...,n— 1,  which  are  mutually  independent. 


where  b  is  a  realization  of  B  =  {Bq,  Hi,  ... ,  Bn-if' ,  p  =  (po,  Pi,  ■  ■  ■ ,  Pn-i)'^  is  a  non-random 
vector  functions  of  t,  7  is  a  non-random  function  of  t,  and  c  is  a  given  constant.  Hence,  it  is 
easy  to  see  that  the  solution  to  the  first  equation  of  (3.16)  with  initial  condition  a;(0)  =  xq 
is  given  by 

x{t;xo,h)  =  —c+{xo  +  c)exp^-J  p{s)ds  +  j  '-y{s)d^  .  (3-17) 

Assume  that  all  the  subsystems  have  the  same  probability  density  function  for  initial  condi¬ 
tion  Xq,  which  is  independent  of  Bj,  j  =  0, 1,  2, . . . ,  n  —  1.  Let  X{t)  =  x{t-,  Xq,  B)  and 


Then  we  have 


y'(t)  =  B  ■  [  p{s)ds  +  [  'y{s)ds. 
Jo  Jo 

X{t)  =  -c+  (Xo  -f  c)exp(F(t)). 


(3.18) 
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Since  Bj  ~  cr‘^),  and  Bj,  j  =  0, 1,  2, . . . ,  n  —  1,  are  mutually  independent,  we  find  that 

{H- p{s)  +  j{s))ds,  pj{s)ds 

Stochastic  formulation  By  the  above  discussions  on  the  probabilistic  formulation  (3.16), 
we  see  that  if  all  the  individuals  in  the  entire  system  have  the  same  fixed  initial  condition  xo, 
then  for  any  fixed  t,  X{t)  has  a  log-normal  distribution  when  c  =  0,  and  a  shifted  log-normal 
distribution  when  c  is  nonzero.  Based  on  this  information,  the  stochastic  model  is  chosen  to 
have  the  form 


dX{t)  =  ^{t){X{t)  +  c)dt  +  ri{t){X{t)  +  c)dW{t),  X(0)  =  Xq,  (3.20) 

where  both  ^  and  rj  are  non-random  functions  of  t,  and  Xq  is  independent  of  W{t).  Let 
h{x,  t)  =  ln(a;  -f  c).  Then  we  have 

1  1 
ht{x,t)  =  0,  h^{x,t)  =  — ^ — ,  h^^{x,t)  =  -- — ^ 

X  +  C  [X  +  C)^ 

Hence,  by  Ito’s  formula,  we  have 


dh{X{t),t)  =  {^{t)  —  ^ri^{t))dt  +  ri{t)dW{t). 
Integrating  both  sides  we  have 

X(t)  =  -c+  (Xo  -h  c)exp(F(t)). 


(3.21) 


where  Y (t)  =  /  (^(s)  —  ■^ri^{s))ds  +  /  ri{s)dW{s).  By  Lemma  3.2,  for  any  fixed  t  we  have 
Jo  Jn 

that 


'0 


r]{s)dW{s)  JY  J  r]^{s)ds^  . 


Hence,  we  find  that  for  any  fixed  t  we  have 

rt 


H(t)  ~  A/"  (^(s)  -  \rf{s))ds,  T]‘^{s)d.s^  . 


(3.22) 


Equivalent  conditions  By  (3.18),  (3.19),  (3.21),  and  (3.22),  we  see  that  if  functions  rj 
and  pj,  and  constants  pj,  aj  and  n  satisfy  the  following  two  equalities 


\p^{s))ds  = 


{p  ■  p{s)  +  ^{s))  ds 


(3.23) 


and 


(3.24) 
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then  probabilistic  formulation  (3.16)  and  stochastic  formulation  (3.20)  yield  processes  that 
are  pointwise  equivalent  in  density. 

Based  on  the  above  equivalent  conditions,  we  next  derive  the  specihc  forms  of  the  func¬ 
tional  parameters  and/or  values  of  the  parameters  of  the  corresponding  pointwise  equivalent 
(stochastic/probabilistic)  formulation  in  terms  of  the  functional  parameters  and/or  param¬ 
eters  of  the  known  (probabilistic/stochastic)  formulation. 


3.2.1  Probabilistic  Formulation  to  Stochastic  Formulation 

Here  we  assume  that  probabilistic  formulation  (3.16)  is  known,  and  we  wish  to  determine 
its  corresponding  pointwise  equivalent  stochastic  formulation.  Differentiating  both  sides  of 
(3.24)  with  respect  to  t  we  have 

n— 1 


=  2  5^a|pj(t)  pj{s)d^  . 


Hence,  if  we  assume  that  Pj(t),  j  =  0, 1,  2, . . . ,  n  —  1  have  the  property  that  the  right-hand- 
side  of  the  above  equation  is  nonnegative  for  any  f  >  0,  then  we  can  always  hnd  p  such  that 
(3.24)  holds,  and  it  is  given  by 


n— 1 


2  W  ( L  Pj{s)ds 

.  j=o 


pit)  = 

By  (3.23)  and  (3.24)  we  hnd  that 

nt  pt  ^  1  /  pt 

/  ^is)ds=  /  (^■p(s) +7(s))ds  + -  (  /  pjis)di 

Jo  Jo  ^  \Jo 

Differentiating  both  sides  of  the  above  equation  with  respect  to  t  yields 

n-l  / 


(3.25) 


^{t)  =  pt  ■  p{t)  +  7(t)  +  J2(r]pjit)  pjis)ds^ 


(3.26) 


Thus,  if  ^  and  p  are  chosen  as  those  in  (3.26)  and  (3.25),  then  both  (3.23)  and  (3.24)  hold. 


3.2.2  Stochastic  Formulation  to  Probabilistic  Formulation 


Now  we  assume  that  stochastic  formulation  (3.20)  is  known,  and  we  wish  to  determine  its 
corresponding  probabilistic  formulation.  It  is  obvious  that  if  we  set 

pt-p(t)  +  7(t)  =^(t)-ip2(t),  (3.27) 

then  (3.23)  holds.  By  (3.27)  and  (3.24)  we  know  that  we  have  lots  of  different  choices  for 
the  probabilistic  formulation.  Here  we  illustrate  two  of  them. 
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Choice  1  One  simple  choice  is  to  choose 

n  =  2,  7  =  0,  /io  =  1,  hi  =  0,  po{t)  =i{t) (3.28) 

Then  by  the  above  equalities,  (3.24)  and  (3.27)  we  obtain  that 

Pi(s)dsj  =  rf{s)ds  -  al  (^(s)  -  \rf{s))d^  ,  (3.29) 

which  implies  that  we  need  to  choose  ag  sufficiently  small  such  that  its  right-hand  side  is 
greater  than  0.  Now  by  (3.29)  we  have 

Tf{s)ds  -  al  (^(s)  -  l7]^(s))ds 

Differentiating  both  sides  of  above  equation  with  respect  to  t  we  obtain  that 

^  ^  p^(s)ds  -  al  (^(s)  -  ^r]‘^{s))ds^ 

Hence,  we  can  just  assign  any  positive  value  for  ai,  and  then  use  the  above  equality  to 
determine  function  pi.  Thus,  one  simple  choice  of  the  corresponding  pointwise  equivalent 
probabilistic  formulation  for  the  stochastic  formulation  (3.20)  is  as  follows 

^  ~  ^0,  bi)  +  c),  bo,  bi  e  M, 

Ho  ~  A/”(l,  (Jg),  Hi  ~  A/”(0,  (Ti), 

where  ag,  ai  and  pi  are  chosen  according  to  the  above  discussions. 


(3.30) 

(3.31) 


then  (3.30)  holds.  Thus,  another  simple  choice  of  the  corresponding  pointwise  equivalent 
probabilistic  formulation  for  the  stochastic  formulation  (3.20)  is  as  follows 

=  {bpo{t)  +  ^{t)-  Wit))  (xit;  b)+c),  beR, 

B  ~  A/'(0, 1), 

where  function  pg  is  dehned  by  (3.31). 
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3.3  System  of  Differential  Equations 


In  this  section,  we  derive  the  conditions  for  the  systems  under  which  the  probabilistic  and 
stochastic  formulations  have  the  same  distribution  (multivariate  normal  in  the  case  the  initial 
condition  is  a  constant  vector)  at  each  time  t. 


Probabilistic  formulation  The  probabilistic  formulation  considered  has  the  following 
form 

^  >4(t)x(t;  b)  +  7(t)  +  7f(t)b, 

B  ~  A/” (/Xg,  Sb), 

where  Al  is  a  non-random  mxm  matrix  functions  of  f,  7  is  a  non-random  m  vector  functions 
of  t,  is  a  non-random  mxn  matrix  functions  of  t,  and  b  is  a  realization  of  B.  The  solution 
to  the  hrst  equation  of  (3.32)  with  initial  conditions  x(0)  =  xq  is  given  by 


x(t;  xo,  b)  =  #(t)xo  -h  $(t) 


r*f 


#  ^{s)j{s)ds+  /  $  ^(s)7-f(s)bds 

lo  Jo 


(3.33) 


Here  $(t)  is  the  solution  of  deterministic  initial  value  problem 

^  =  #(0)=I„,  (3.34) 

where  Xm  is  m  x  m  identity  matrix.  Assume  that  all  the  subsystems  have  the  same  proba¬ 
bility  density  function  for  initial  condition  Xq,  which  is  independent  of  B,  and  let  X(t)  = 
x(f;Xo,  B)  and 


Y(t)  =  <!>(«) 


^{s)^{s)ds  +  [  ^{s)H{s)'Bds 

Jo 


Then  we  have  that 


X(t)  =  $(t)Xo  +  Y(t). 

Note  that  B  ~  Sb).  Hence,  we  hnd  that  for  any  hxed  t  we  have 

rt 


Y{t)r^Af~  $(t)  /  ^-\s)  +  n{s)ti^)ds,^{t)np{t){^{t)y 


where 


Up{t)  = 


^~\s)n{s)ds 


L^o 


r-i 

^~\s)n{s)ds  . 


L^o 


(3.35) 


(3.36) 
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Stochastic  formulation  By  the  above  discussions  on  the  probabilistic  formulation  (3.32), 
we  see  that  if  all  the  individuals  in  the  entire  system  have  the  same  hxed  initial  condition 
xq,  then  X(t)  is  multivariate  normal  distributed  for  any  hxed  time  t.  Based  on  this  piece  of 
information,  the  stochastic  model  is  chosen  to  have  the  following  form: 

dX(t)  =  [A{t)X{t)  +  ^{t)]dt  +  J^{t)dW{t),  X(0)  =  Xo,  (3.37) 


where  ^  is  a  non-random  m  vector  functions  of  t,  JF  is  a  non-random  m  x  /  matrix  function  of 
t,  and  W(t)  =  (hhi(t),  13^2 (t), . . . ,  Wi{t)Y'  is  a  /-vector  standard  Wiener  process  independent 
of  initial  vector  Xq.  Then  the  solution  to  (3.37)  is  given  by  (e.g.,  see  [26,  Section  5.6]) 


X(/)  =  $(/)Xo  +  Y(/), 


where  Y (/)  is  dehned  by 


Y{t)  =  ^{t)  [  ^-\s)^{s)ds+  [  ^~\s)J^{s)dW{s 

Uo  Jo 

In  addition,  for  any  hxed  t  we  have  (e.g.,  see  [26,  Section  5.6]) 

Y(«)  ~  M  (^^t)  #(()n.(«)(#(()f])  , 


where 


n,(t)=  [  ^-\s)J^{s){^-\s)J^{s)fds. 

Jo 


(3.38) 


(3.39) 


Equivalent  conditions  By  (3.35),  (3.36),  (3.38)  and  (3.39),  we  know  that  if  functions 
JF,  7,  H,  /Xg  and  Sb  satisfy  the  following  two  equalities 


#  {s){'y{s)  +n{s)fie)ds  =  /  #  {s)^{s)ds, 


(3.40) 


and  np(/)  =  n5(/),  that  is. 


^-\s)nis)ds 


->T 


^-\s)n{s)ds 


=  /  ds, 


(3.41) 

then  probabilistic  formulation  (3.32)  and  stochastic  formulation  (3.37)  are  pointwise  equiv¬ 
alent  in  density. 


Next  we  will  derive  the  specihc  forms  of  the  functional  parameters  and/or  values  of  the  pa¬ 
rameters  of  the  corresponding  pointwise  equivalent  (stochastic/probabilistic)  formulation  in 
terms  of  the  functional  parameters  and/or  parameters  of  the  known  (probabilistic/stochastic) 
formulation  based  on  the  equivalent  conditions  (3.40)  and  (3.41). 
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3.3.1  Probabilistic  Formulation  to  Stochastic  Formulation 


Here  we  assume  that  probabilistic  formulation  (3.32)  is  known,  and  we  want  to  determine 
its  corresponding  pointwise  stochastic  formulation.  By  (3.40),  it  is  obvious  that  if  function 
^  is  chosen  to  be 

then  (3.40)  holds.  Differentiating  both  sides  of  (3.41)  with  respect  to  t  we  hnd  that 

T(t)T'^(t)  =  <I.(i)n,(i)(<l>(«))^.  (3.42) 

Hence,  if  we  assume  lip  have  the  property  that  #(t)np(f)($(t))^  is  a  positive-semidehnite 
matrix  for  any  f  >  0,  then  we  can  always  hnd  T  such  that  (3.42)  holds. 


3.3.2  Stochastic  Formulation  to  Probabilistic  Formulation 

Now  we  assume  that  stochastic  formulation  (3.37)  is  known,  and  we  want  to  determine  its 
corresponding  pointwise  equivalent  probabilistic  formulation.  It  is  obvious  if  7,  and 
are  chosen  such  that 

7(t)  +  7f(t)ArB  = 

then  (3.40)  holds.  Note  that  Ti  also  needs  to  be  properly  chosen  to  satisfy  (3.41).  Hence, 
for  simplicity,  we  choose 

7(t)=|(t),  /Xb  =  0.  (3.43) 

Note  that  for  any  f  >  0,  n5(t)  is  a  positive-semidehnite  matrix.  Hence,  there  exists  a  mx  m 
matrix  function  A{t)  such  that 


n.(t)  =  A(t)(A(t))^. 


Thus,  if  we  choose 

n  =  m,  H{t)  =  (3.44) 

then  (3.41)  holds.  Therefore,  one  of  the  simple  choices  of  the  corresponding  pointwise  equiv¬ 
alent  probabilistic  formulation  for  the  stochastic  formulation  (3.37)  is  as  follows 


=  A(t)x(f;  b)  +  |(t)  +  $(t)  A(f)b,  b  G  W^, 

B  ~  Af{0,ljn)- 


4  Computing  Covariances 


As  we  have  already  stated,  the  procedures  above  lead  to  pointwise  equivalent  in  density 
systems  that  are  not  described  by  the  same  stochastic  processes.  We  revisit  the  examples 
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of  (2.1)  to  illustrate  this  by  comparing  covariances  of  the  two  processes.  We  proceed  to  use 
the  Lemmas  3.1  and  3.2  to  hnd  the  covariance  function  of  the  stochastic  processes 
in  the  probabilistic  formulation  (2.4)  and  YsRuit)  in  the  stochastic  formulations  (2.5). 

Probabilistic  formulation:  Note  that  for  this  case  Yppp  is  given  by 

YpRuit)  =  exp{Bt  -  where  B  ~  M{bo,  Uq). 

Hence,  by  Lemma  3.1  we  hnd  immediately 

E(ypijD(t))  =  exp(6ot).  (4.1) 

Then  using  Lemma  3.1  and  (4.1)  we  hnd  the  covariance  function  for  the  process  {H(t)}  = 


{YpRB(t)}  given  by 

Cov(Y(t),Y(s)) 

=  E(Y(t)Y(s))  -  E(Y(t))E(Y(s)) 

=  E  {exp  (B(t  +  s)  -  |(To(t^  +  s^))  }  -  exp(6o(t  +  s)) 

=  exp  (bo(t  +  s)  -  |(To(t^  +  s^)  +  |(To(t  +  s)^)  -  exp(&o(t  +  s)) 
=  exp  (6o(t  +  s)  +  st(To)  -  exp(6o(t  +  s)) 

=  exp(6o(t  +  s))  [exp  (stag)  -  l]  . 

Stochastic  formulation:  In  this  case  we  have 


YsRMit)  =  exp  (^{bot  -  |cTot^)  +  j  V^dWir)^ 


Let  Q{t)  =  (Jo  V^dW{T).  Then  by  Lemma  3.2,  we  have  that  {(^(t)}  is  a  Gaussian 

Jo 

process  with  zero  mean  and  covariance  function  given  by 

Cov{Q{t),  Q(s))  =  (To  min{t^,  s^}.  (4.2) 

Using  Lemma  3.1  and  (4.2)  we  hnd  that 


E(YsRM(t))  =  exp(bot). 


(4.3) 


By  Lemma  3.2  we  know  that  {Q(t)}  is  a  Gaussian  process.  Hence,  Q(t)+Q(s)  has  a  Gaussian 
distribution  with  zero  mean  and  variance  dehned  by 


Var(Q(t)  +  Q(s)) 


Var(Q(t))  +  Var((5(s))  +  2Cov(Q(t),  Q(s)) 
(Tg  (t^  +  +  2  min{t^,  s^})  . 


(4.4) 
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Now  we  use  Lemma  3.1,  along  with  equations  (4.3)  and  (4.4)  to  find  the  covariance  function 
of  {Yit)}  =  {YsRuit)}. 

Cov{Y{t),Y{s)) 

=  EiYit)Yis))-EiYit))EiYis)) 

=  E  {exp  (bo{t  +  s)  -  +  s^)  +  Q{t)  +  Q(s)) }  -  exp(6o(t  +  s)) 

=  exp  (bo(t  +  s)  -  |no(t^  +  s^)  +  +  2  min{t^,  s^}))  -  exp(6o(t  +  s)) 

=  exp  (bo(t  +  s)  +  (Tq  min{t^,  s^})  -  exp(6o(t  +  s)) 

=  exp(6o(t  +  s))  [exp  min{f^,  s^})  -  l]  . 

Thus  we  see  that  {YpR^it)}  and  {YsRMit)}  have  different  covariance  functions  and  hence 
are  not  the  same  stochastic  process. 


5  Equivalence  between  Probabilistic  and  Stochastic  For¬ 
mulations  with  Nonlinear  Dynamics 


In  this  section  we  turn  to  nonlinear  stochastic  differential  equations  which  can  be  shown 
equivalent  to  a  PRD  formulation  with  nonlinear  dynamics.  In  summary  of  our  results  to  date, 
based  on  the  discussions  in  Sections  3.1  and  3.2,  we  see  that  we  can  find  the  corresponding 
pointwise  equivalent  probabilistic  formulation  for  two  types  of  scalar  stochastic  differential 
equations: 

dX{t)  =  [a{t)X{t)  +  ^{t)]dt  +  rj{t)dW{t), 

and 

dX{t)  =  i{t){X{t)  +  c)dt  +  ri{t){X{t)  +  c)dW{t), 

where  t],  and  a  are  all  deterministic  functions  of  t,  and  c  is  a  given  constant.  In  addition, 
we  can  find  the  corresponding  pointwise  equivalent  probabilistic  formulation  for  the  following 
linear  system  of  stochastic  differential  equations  (see  Section  3.3) 

dX{t)  =  [A{t)X{t)  +  ^{t)]dt  +  X{t)dW{t), 

where  ^  is  a  non-random  mxm  matrix  functions  of  f,  ^  is  a  non-random  m  vector  functions 
of  t,  JF  is  a  non-random  m  x  I  matrix  function  of  t,  and  W(t)  =  W2{t), . . . ,  Wi{t))'^  is 

a  /-vector  standard  Wiener  process.  Hence,  if  a  nonlinear  stochastic  differential  equation  (or 
system  of  stochastic  differential  equations)  can  be  reduced  to  one  of  the  above  forms  by  some 
invertible  transformation,  then  one  can  find  its  corresponding  probabilistic  formulation.  The 
same  thing  is  true  for  the  probabilistic  formulation. 
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5.1  Reducible  Nonlinear  Stochastic  Differential  Equations 

First  we  will  consider  several  special  cases  of  nonlinear  stochastic  differential  equations  that 
can  be  reduced  to  linear  stochastic  differential  equations  after  some  invertible  transformation. 
Proofs  are  given  in  [24,  Section  4.1], 


Theorem  5.1  Consider  the  stochastic  differential  equation 

dX{t)  =  g{X{t)C)dt  +  a{X{t)C)dW{t), 
where  g  and  a  are  non-random  functions  of  x  and  t.  If  the  equality 


A 

dx 


a(x,  t) 


1  da  d{g\  1 3“^  a 


=  0 


holds,  then  the  nonzero  deterministic  function  a(t)  can  be  determined  from 

1  da ,  d  f  g\  ,  1  dd‘a , 


a'{f)  =  a{t)a{x,  t) 


■  /  N  f  Q\  ,  N  t  o  a ,  , 


la^(x,  f)  dt 

and  some  smooth  invertible  function  h{x,  f)  can  be  computed  from 

dh ,  ,  dif) 

dx  a{x,  t) 

Moreover,  (5.1)  can  be  reduced  to  the  linear  stochastic  differential  equation 

dZ{f)  =  g{t)dt  +  d{t)dW  it), 

where  Z(t)  =  h{X(t),t)  and  the  deterministic  function  gif)  can  be  computed  from 

-/  N  dh ,  ,  dh ,  ,  ,  ,  13‘^h,  ,9,  , 

®  ^  Ht  ^  dx^^'  *  ®  ^  2 *)■ 


(5.1) 


Theorem  5.2  The  autonomous  stochastic  differential  equation 

dXit)  =  g{X{t))dt  +  a{X{t))dW{t) 
can  be  reduced  to  the  linear  stochastic  differential  equation 

dZit)  =  (Ao  +  \iZ{t))dt  +  (z/o  +  viZ{t))dW{t) 


if  and  only  if 


fjfx)  =  0  or 


[af^ 

if}' 


(x)  =  0. 


(5.2) 


Here  g  and  a  are  non-random  functions  of  x,  Xq,  \i,  uq  and  z/i  are  some  constants,  ipi^x)  = 

g  1 

— - a'ix),  and  Zif)  =  hiXif)),  where  h  is  some  invertible  transformation. 

a{x)  2 
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Ip' 


If  vi  7^  0,  then  we  ean 


If  the  latter  part  of  (5.2)  is  satisfied,  then  we  set  ui  = 
choose 


where  c  is  some  constant. 


h{x)  =  cexp  (  vi 


O'  r 


If  ui  =  0,  then  we  can  choose 


1 


h{x)  =  uq 


1 

(j{t) 


dr  +  c. 


Remark  5.3  For  a  general  nonlinear  system  of  stochastic  differential  equations,  it  is  difficult 
to  obtain  the  explicit  form  for  the  invertible  transformation  (as  it  strongly  depends  on  the 
specific  form  of  the  system).  Hence,  we  do  not  pursue  this  effort  in  this  paper. 


5.2  Examples 

We  next  use  several  examples  to  illustrate  this  transformation  method  to  hud  the  correspond¬ 
ing  equivalent  probabilistic/stochastic  formulations  for  stochastic/probabilistic  formulations. 


Example  5.1  (Exponential  Modulated  Growth  in  the  Drift) 

In  this  example,  we  use  the  transformation  method  to  hud  the  equivalent  probabilistic  for¬ 
mulation  for  the  following  nonlinear  stochastic  differential  equation 


dX  = 


1  -  2 


dt  -|-  exp{—X)dW. 


Note  here  that  g{x)  ~  2  ^  ~  exp(X).  Then  by  Ito’s  formula  we  hnd  that 


dZ  =  exp(X) 


1  -  2 


■1 


dt  -|-  exp(— X)(iIT  J  +  2  exp(— 2X)(it 


=  Zdt  +  dW. 


By  the  discussions  in  Section  3.1.2  we  hnd  the 
above  linear  stochastic  differential  equation  is 


dzit]  b) 
dt 


=  z{t;b)  +bp{t), 


where  pit)  = 


exp(2t) 


a/2  [exp  (2t)  —  1] 


exp(2t) 


dx 

dt 


^{z  +  bp{t))  = 


equivalent  probabilistic  formulation  for  the 

b  eR;  B  A/'(0, 1), 

Let  X  =  ln(;2).  Then  we  have 
1  -|-  bp(t)  exp(— x). 
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Thus,  the  following  two  formulations 


dX 


1 


1 

2 


exp(— 2X) 


dt  +  exp{—X)dW 


and 


dx(t;b)  ^  exp(2t) 

dt  A/2|exp(2ty^^TJ 

B  ~  A/'(0, 1) 

are  pointwise  equivalent  in  density. 


exp(2f)  —  1 


exp(-a;(f;6)), 


6  G  M, 


Example  5.2  (Logistics  Growth  Dynamics  with  Uncertainty) 

We  begin  with  logistic  growth  in  the  probabilistic  formulation  and  derive  its  equivalent 
stochastic  formulation.  Consider  the  deterministic  logistic  equation 

^  =  hx(l-^  ,  a;(0)  =  xq,  (5.3) 

where  b  is  some  constant  representing  the  intrinsic  growth  rate,  and  k  is  a  given  constant 
representing  the  carrying  capacity.  Let  =  — .  Then  it  is  easy  to  find  that 

X 


dz 

dt 


By  the  discussions  in  Section  3.2.1,  we  know  that  for  probabilistic  formulation 

=  -b  (^zit]  b)-^  ,  beR;  B  A/'(/io,  (Tq), 


its  equivalent  form  of  stochastic  formulation  is  given  by 

dZ{t)  =  (— po  +  CTQt)  - ^  dt  +  v^uo  ^Z 


K 


dW{t). 


Let  X{t)  = 


Z{t) 


Then  by  Ito’s  formula  we  find  that 


dX{t)  = 


Z^it)  [ 
1 


(-fto  +  crlt)  ( Z  -  dt  +  (z--]  dW {t) 


(  Z - 

V  « 


dt 


Z%t)  [ 

=  X  alt)  +  2tal 


dt  —  V^aoX  - ^  dW(t). 


25 


Thus,  the  probabilistic  formulation 

=  bx(t;  b)  ,  6  G  M;  S  ~  Afi.Ho,  a^)  (5.4) 

and  the  stochastic  formulation 

^  dW{t) 

(5.5) 

are  pointwise  equivalent  in  density.  Figure  1  depicts  the  probability  density  function  p{x,t) 
at  different  times  t  for  the  probabilistic  formulation  (5.4)  and  the  stochastic  formulation 
(5.5)  with  K  =  100,  xq  =  10,  po  =  1  and  ag  =  0.1,  where  p{x,t)  is  obtained  by  simulating 
10^  sample  paths  for  each  formulation.  Here  for  the  probabilistic  formulation  (5.4),  we 


dX{t)  =  X 


(po  -  o-Qt)  (  1  -  ^  )  +  2tal 


X 

1 - 

K 


dt  —  V^agX  (1 - 

V 
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Figure  1:  Probability  density  function  p{x,t)  are  obtained  by  simulating  10®  sample  paths 
for  probabilistic  formulation  (5.4)  and  stochastic  formulation  (5.5)  at  t  =  1,  2,  3  and  4,  where 
At  =  0.004  is  used  in  (5.6),  and  T  =  4. 

analytically  solve  each  deterministic  differential  equation  (5.3)  with  b  being  a  realization  of 
B,  which  follows  a  normal  distribution  A/'(l,  0.01).  This  solution  is  given  by 

^  Xo^exp(5t) 

’  K  +  a;o(exp(6t)  —  1) 
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Thus  implementation  of  the  probabilistic  formulation  (PRD)  is  extremely  rapid  (essentially 
a  function  evaluation).  This  is  not  true  for  the  SRM  and  hence  the  PRD  formulation 
requires  much  less  (orders  of  magnitude)  implementation  time  than  does  the  SRM.  Even  in 
examples  where  both  methods  require  numerical  integration,  the  PRD  is  highly  preferred 
to  the  SRM  in  implementation  of  inverse  problems  (e.g.,  see  [11,  14,  20]).  We  use  an  Euler 
explicit  method  (crude  but  sufficient  for  our  purposes  here  since  our  main  objective  is  to 
demonstrate  our  theoretical  equivalence  results  with  a  numerical  example)  to  numerically 
approximate  the  sample  paths  for  the  stochastic  differential  equation  (5.5).  Let  T  denote 
the  hnal  time,  and  m  be  the  number  of  mesh  points  interval.  Then  the  mesh  time  points 
are  given  by  tk  =  kAt,  k  =  0, 1,  2, . . . ,  m,  where  At  =  T /m.  Denote  by  the  numerical 
solution  for  X{tk)',  then  we  have  the  following  numerical  scheme: 


^  x’^  +  X>^  (/io-a24)(l 


X^ 


X^ 


K 


+  2tfc(TQ  I  1 


X 


fc\  2' 


K 


At 


(5.6) 


-v/24aoX^  (1 - )  £k,  k  =  1,2, . . .  ,m  -  1. 

n 


where  Sk  is  a  random  variable  following  a  normal  distribution  A/'(0,  At).  From  Figure  1  we 
see  that  we  obtain  the  same  probability  density  function  for  the  probabilistic  formulation 
(5.4)  and  the  stochastic  formulation  (5.5),  which  nicely  illustrates  our  earlier  theoretical 
results. 


Example  5.3  (Gompertz  Growth  in  the  Drift) 

In  this  example,  we  consider  the  following  nonlinear  stochastic  differential  equation 

dX{t)  =  [ao{t)  -  ai(t)  ln(X(t))]  X{t)dt  +  y^2d^X{t)dW{t),  (5.7) 


where  Oq,  Oi  and  do  are  some  deterministic  functions  of  t,  and  Oi  and  do  ai’e  assumed 
to  be  positive.  This  equation  is  a  stochastic  version  of  the  generalized  Gompertz  model 
X  =  (ao(t)  —  ai(t)  lna;)a;,  which  has  been  extensively  used  in  biological  and  medical  research 
to  describe  population  dynamics  such  as  tumor  growth  in  humans  and  animals  either  with 
or  without  treatment  (e.g.,  [1,  23]  and  the  references  therein)  as  well  in  cell  proliferations 
models  [18].  Next  we  use  the  transformation  method  to  hnd  an  equivalent  probabilistic 
formulation  for  (5.7).  Let  Z  =  ln(X).  Then  by  Ito’s  formula  and  (5.7)  we  have 


dZ{t) 


|[ao(t)  -ai{t)  ln{X {t))]  X {t)dt  +  ^/2dJ^X{t)dW{t)^ 
[2do(t)X2(t)]  dt 

[-ai{t)Z{t)  +  (ao(t)  -  do(f))]  dt  +  ^/2do{t)dW (t) . 


By  the  discussions  in  Section  3.1.2  we  hnd  that  a  pointwise  equivalent  probabilistic  formu¬ 
lation  for  the  above  linear  stochastic  differential  equation  is  given  by 


=  -ai{t)z{t;  b)  +  ao{t)  -  do{t)  +  bp{t),  b  eR, 
B  ~  A/'(0, 1). 
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Here  6  is  a  realization  of  B,  and  p  is  given  by 

d 


pit)  =  ^ 


where 

(pit)  =  I  2(io(-s)exp  (  —2  j  ai(r)(ir  )  ds. 

Jo 

Let  X  =  exp(;2).  Then  by  (5.8)  we  hnd 
dxit;  b) 


dt 


=  exp{z{t]  b))  [-ai{t)z{t]  b)  +  ao(t)  -  doit)  +  &p(t)] 


=  a;(t;  b)  [-ai(t)  ln(a;(t;  b))  +  ao(t)  -  doit)  +  &p(t)] . 
Thus,  the  following  two  formulations 


dXit)  =  [ao(t)  -  ai(f)  ln(X(f))]  Xit)dt  +  ^2doit)Xit)dWit) 

and 

^^^dt  ^  \-~^^it)  ^))  +  ®o(^)  —  doit)  +  bpit)] ,  &  G  M 

B  ~  A/'(0, 1) 

are  pointwise  equivalent  in  density. 


6  Concluding  Remarks 


In  summary,  we  have  derived  several  classes  of  examples  with  affine  dynamics  for  which 
we  can  establish  pointwise  equivalence  in  density  for  the  corresponding  probabilistic  and 
stochastic  formulations.  We  then  argue  that  a  large  class  of  nonlinear  SDE  can  be  reduced 
by  invertible  transformation  to  one  of  the  affine  cases.  We  presented  several  examples  of  non¬ 
linear  SDE  arising  frequently  in  applications  that  can  be  transformed  to  the  more  readily 
computed  probabilistic  formulation.  It  has  been  well  documented  (e.g.,  see  [20,  25]  and  ref¬ 
erences  therein)  that  difficulties  arise  in  numerically  solving  Fokker-Planck  equations  such  as 
(1.1)  when  the  drift  g  dominates  the  diffusion  This  motivated  our  efforts  and  the  results 
here  lead  to  alternative  methods  that  can  be  fast  and  efficient  in  numerically  solving  (1.1) 
by  employing  its  pointwise  equivalent  in  density  probabilistic  formulation.  Future  efforts  in¬ 
clude  investigation  of  the  approximate  equivalent  probabilistic  formulations  for  non-reducible 
nonlinear  stochastic  differential  equations  by  the  linearization  method  as  well  investigation 
of  applications  of  our  methodology  to  the  control  of  systems  with  uncertainty. 
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